for(i in 1:17){
  
VI <- dependent_vars[i]
mod <- readRDS(paste0('./models/',VI,'_poly2.rds'))

print(VI)
print(pp_check(mod))
print(mcmc_plot(mod,type = "rhat"))
}
## [1] "CCI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "CIgr"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "CIre"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "NDVI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "NRVIre"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "SR"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "ARI1"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "ARI2"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "CRI1"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "CRI2"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "PRI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "MSI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "NDWI1"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "NDWI2"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "SRWI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "WI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## [1] "WI_NDVI"
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

for(i in 1:17){
  mod <- readRDS(paste0('./models/',dependent_vars[i],'_poly2.rds'))
 
em_week <- emmeans(
  mod, 
  specs = pairwise ~ treatment | days_scaled,
  at = list(days_scaled = seq(min(data$days_scaled), 
                              max(data$days_scaled), 
                              length = 8))  # Evaluates at original weeks
)

# Plot treatment differences
p1 <- plot(em_week$contrasts) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Drought vs. Control: Weekly Contrasts",
       y = "VI Difference (Drought - Control)")

em_clone <- emmeans(
  mod, 
  specs = pairwise ~ treatment | days_scaled| clone,
  at = list(days_scaled = seq(min(data$days_scaled), 
                              max(data$days_scaled), 
                              length = 8))  # Evaluates at original weeks
)

contrasts_df <- as.data.frame(em_clone$contrasts)
p2 <- ggplot(contrasts_df, aes(x = days_scaled, y = estimate, color = contrast)) +
  geom_point() +
  geom_line() +
  geom_errorbar(aes(ymin = lower.HPD, ymax = upper.HPD), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Drought vs. Control: Weekly Contrasts by Clone",
       y = "VI Difference (Drought - Control)",
       x = "Days Scaled") +
  facet_wrap(~ clone) +  # Facet by clone
  theme_minimal() +
  #annotate("text", x=1, y=1, label= annotate("text", label= dependent_vars[i]))+
  theme(legend.position = "bottom")
 print(dependent_vars[i])
print(p1)
print(p2)
}
## NOTE: Results may be misleading due to involvement in interactions
## [1] "CCI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "CIgr"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "CIre"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "NDVI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "NRVIre"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "SR"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "ARI1"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "ARI2"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "CRI1"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "CRI2"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "PRI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "MSI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "NDWI1"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "NDWI2"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "SRWI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "WI"

## NOTE: Results may be misleading due to involvement in interactions
## [1] "WI_NDVI"